import pandas as pd
pd.set_option('display.expand_frame_repr', False)

import matplotlib
import matplotlib.pyplot as plt
font = {'size' : 20}
matplotlib.rc('font', **font)

from pylab import rcParams
rcParams["figure.figsize"] = 30,16

import seaborn as sns

from statsmodels.tsa.statespace.sarimax import SARIMAX
from pmdarima.arima import auto_arima

import datetime as dt
from datetime import timedelta, date

from sklearn.preprocessing import MinMaxScaler
import joblib

import warnings
warnings.filterwarnings("ignore")

import sys
sys.path.insert(0, "../")
import functions

#Unskalierte Daten für Analysen laden
df_unscaled = pd.read_csv("../3-Data Preparation/arima_data_unscaled.csv", index_col=0, parse_dates=True)
df_unscaled.index.freq = "D"

#Wahre Werte der Testdaten
y_true = df_unscaled["verbrauch"]["2021-01-01":]

#Kalender laden
df_calendar = pd.read_csv("../2-Data Understanding/Datenbeschaffung/kalender.csv", index_col=0, parse_dates=True)

#Skalierte Daten für Modellierung laden
df_scaled = pd.read_csv("../3-Data Preparation/arima_data_scaled.csv", index_col=0, parse_dates=True)
df_scaled.index.freq = "D"

#Aufteilung in endogene und exogene Daten
exog_train = df_scaled[["arbeitstag", "temperatur"]][:"2020-12-31"]
exog_test = df_scaled[["arbeitstag", "temperatur"]]["2021-01-01":]

exog = exog_train.append(exog_test)

endog_train = df_scaled["verbrauch"][:"2020-12-31"]
endog_test = df_scaled["verbrauch"]["2021-01-01":]

endog = endog_train.append(endog_test)

Analyse der besonders großen Abweichungen

Das ARIMA(2,0,2)(2,0,2)7 [Arbeitstag, Temperatur] hat große Probleme an und nach Feiertagen. Ziel ist es, diese Abweichungen durch weitere Merkmale zu verringern.

ARIMA(2,0,2)(2,0,2)7 [Arbeitstag, Temperatur]

#Nicht-saisonale und saisonale Ordnung festlegen
order = (2, 0, 2) #p, d, q
seasonal_order = (2, 0, 2, 7) #P, D, Q, m

#Modell mit Trainingsdaten erstellen
train_model_config = SARIMAX(endog=endog_train, exog=exog_train, order=order, seasonal_order=seasonal_order)
train_model = train_model_config.fit()

#Modell mit allen Daten erstllen und Konfiguration/Koeffizienten von ersten Modell übernehmen (kein neues Training)
model_config = SARIMAX(endog=endog, exog=exog, order=order, seasonal_order=seasonal_order)
model = model_config.filter(train_model.params)

#Modell ausgeben
print(model.summary())
print()

                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                         verbrauch   No. Observations:                 2557
Model:             SARIMAX(2, 0, 2)x(2, 0, 2, 7)   Log Likelihood                4556.901
Date:                           Thu, 03 Feb 2022   AIC                          -9091.801
Time:                                   22:26:17   BIC                          -9027.489
Sample:                               01-01-2015   HQIC                         -9068.479
                                    - 12-31-2021                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
arbeitstag     0.2951      0.002    151.520      0.000       0.291       0.299
temperatur    -0.0734      0.014     -5.418      0.000      -0.100      -0.047
ar.L1          1.1178      0.102     10.921      0.000       0.917       1.318
ar.L2         -0.1766      0.091     -1.950      0.051      -0.354       0.001
ma.L1         -0.2569      0.103     -2.486      0.013      -0.459      -0.054
ma.L2         -0.1455      0.015     -9.447      0.000      -0.176      -0.115
ar.S.L7        0.1746      0.638      0.274      0.784      -1.076       1.425
ar.S.L14       0.8174      0.635      1.287      0.198      -0.428       2.062
ma.S.L7       -0.0991      0.631     -0.157      0.875      -1.335       1.137
ma.S.L14      -0.7655      0.581     -1.317      0.188      -1.905       0.374
sigma2         0.0017   2.85e-05     59.953      0.000       0.002       0.002
===================================================================================
Ljung-Box (L1) (Q):                   0.03   Jarque-Bera (JB):              3671.53
Prob(Q):                              0.85   Prob(JB):                         0.00
Heteroskedasticity (H):               0.91   Skew:                            -0.16
Prob(H) (two-sided):                  0.17   Kurtosis:                         8.86
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

#Vorhersage erzeugen
scaled_preds = model.predict(start="2021-01-01", end="2021-12-31", dynamic=False).values

#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics_arima(y_true, scaled_preds, residuals=True)

      Vorhersage
R2          0.96
MAE       3613.0
MSE   26901261.0
RMSE      5187.0
                
MAPE       2.1 %

Große Abweichungen an und um Feiertage

An Feiertagen ist der Stromverbrauch außergewöhnlich niedrig und die wöchentliche Saisonalität wird unterbrochen. Dennoch handelt es sich nicht um Ausreißer im klassischen (rein mathematischen) Sinne, da die Abweichung einer gewissen Systematik folgen. Außerdem hat das Modell nicht nur Probleme an Feiertagen selbst, sondern auch unmittelbar danach. Wenn die außergewöhnlich niedrigen Stromverbräuche als Eingangsvariablen für die Regressionsgleichung folgender Tage verwendet werden, verzerrt sich auch die Vorhersage für folgende Tage. Dadurch kommt es zu starken Abweichungen an einem, zwei, sieben und vierzehn Tagen nach einem Feiertag.

An Feiertagen

#Feiertage einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date, ymin=0, ymax=1, color="orange", linestyle="-.")

#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics_arima(y_true, scaled_preds, False, False, True, False)

Am ersten Tag nach Feiertagen

#Ersten Tag nach Feiertagen einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date + timedelta(days=1), ymin=0, ymax=1, color="orange", linestyle="-.")
    
#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics_arima(y_true, scaled_preds, False, False, True, False)

Am zweiten Tag nach Feiertagen

#Zweiten Tag nach Feiertagen einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date + timedelta(days=2), ymin=0, ymax=1, color="orange", linestyle="-.")
    
#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics_arima(y_true, scaled_preds, False, False, True, False)

Am siebten Tag nach Feiertagen

#Woche nach Feiertagen einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date + timedelta(days=7), ymin=0, ymax=1, color="orange", linestyle="-.")
    
#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics_arima(y_true, scaled_preds, False, False, True, False)

Am vierzehnten Tag nach Feiertagen

#Zwei Wochen nach Feiertagen einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date + timedelta(days=14), ymin=0, ymax=1, color="orange", linestyle="-.")
    
#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics_arima(y_true, scaled_preds, False, False, True, False)

Verbesserung durch Merkmal für Feiertage

Es soll im Folgenden versucht werden, diese Abweichungen durch Ausgleichskoeffizienten zu kompensieren.

Verbesserung durch Merkmal für Feiertage

Dafür wird zunächst ein Indikator für Feiertage eingefügt. Dadurch lässt sich allerdings keine Verbesserung erzielen. Die Abweichung an und um die Feiertage herum wird dadurch allerdings kaum reduziert. Insgesamt verbessert sich der MAPE nicht. Allerdings hat sich das Modell sehr viel besser an die Trainingsdaten angepasst, die Log-Likelihood steigt auf 4.644 and und das AIC fällt auf -9.264. Es kann hier also von einer leichten Überanpassung ausgegangen werden. Der MAPE verbessert sich sehr wahrscheinlich deshalb nicht, weil das Modell Feiertage selbst bisher relativ gut erkennen konnte. Das Problem ergibt sich eher aus den Feiertagen als Eingabeparameter für die Vorhersage weiterer Tage.

#Unskalierte Daten für Analysen laden
df_unscaled = pd.read_csv("../3-Data Preparation/arima_data_unscaled.csv", index_col=0, parse_dates=True)
df_unscaled.index.freq = "D"

#Wahre Werte der Testdaten
y_true = df_unscaled["verbrauch"]["2021-01-01":]

#Skalierte Daten für Modellierung laden
df_scaled = pd.read_csv("../3-Data Preparation/arima_data_scaled.csv", index_col=0, parse_dates=True)
df_scaled.index.freq = "D"

#Exogene Daten um Feiertag erweitern
df_scaled["feiertag"] = df_calendar["feiertag"].notna().astype(int)

#Aufteilung in endogene und exogene Daten
exog_train = df_scaled[["arbeitstag", "temperatur","feiertag"]][:"2020-12-31"]
exog_test = df_scaled[["arbeitstag", "temperatur","feiertag"]]["2021-01-01":]

exog = exog_train.append(exog_test)

endog_train = df_scaled["verbrauch"][:"2020-12-31"]
endog_test = df_scaled["verbrauch"]["2021-01-01":]

endog = endog_train.append(endog_test)

#Nicht-saisonale und saisonale Ordnung festlegen
order = (2, 0, 2) #p, d, q
seasonal_order = (2, 0, 2, 7) #P, D, Q, m

#Modell mit Trainingsdaten erstellen
train_model_config = SARIMAX(endog=endog_train, exog=exog_train, order=order, seasonal_order=seasonal_order)
train_model = train_model_config.fit()

#Modell mit allen Daten erstllen und Konfiguration/Koeffizienten von ersten Modell übernehmen (kein neues Training)
model_config = SARIMAX(endog=endog, exog=exog, order=order, seasonal_order=seasonal_order)
model = model_config.filter(train_model.params)

#Modell ausgeben
print(model.summary())
print()

                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                         verbrauch   No. Observations:                 2557
Model:             SARIMAX(2, 0, 2)x(2, 0, 2, 7)   Log Likelihood                4643.945
Date:                           Thu, 03 Feb 2022   AIC                          -9263.891
Time:                                   22:26:29   BIC                          -9193.732
Sample:                               01-01-2015   HQIC                         -9238.449
                                    - 12-31-2021                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
arbeitstag     0.2395      0.003     83.391      0.000       0.234       0.245
temperatur    -0.0517      0.012     -4.173      0.000      -0.076      -0.027
feiertag      -0.0752      0.003    -23.864      0.000      -0.081      -0.069
ar.L1          1.0159      0.244      4.167      0.000       0.538       1.494
ar.L2         -0.0759      0.224     -0.339      0.735      -0.515       0.364
ma.L1         -0.2257      0.244     -0.923      0.356      -0.705       0.253
ma.L2         -0.0539      0.034     -1.571      0.116      -0.121       0.013
ar.S.L7        0.1224      0.059      2.081      0.037       0.007       0.238
ar.S.L14       0.8729      0.059     14.899      0.000       0.758       0.988
ma.S.L7       -0.0130      0.052     -0.248      0.804      -0.116       0.090
ma.S.L14      -0.8363      0.046    -18.136      0.000      -0.927      -0.746
sigma2         0.0015    2.4e-05     64.354      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                   2.95   Jarque-Bera (JB):              3928.90
Prob(Q):                              0.09   Prob(JB):                         0.00
Heteroskedasticity (H):               0.92   Skew:                            -0.42
Prob(H) (two-sided):                  0.24   Kurtosis:                         9.01
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

#Feiertage einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date, ymin=0, ymax=1, color="orange", linestyle="-.")

#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics_arima(y_true, scaled_preds, True, False, True, False)

      Vorhersage
R2          0.96
MAE       3613.0
MSE   26901261.0
RMSE      5187.0
                
MAPE       2.1 %

Verbesserung durch Merkmal für vergangene Feiertage (1)

Es wird nun der erste Ausgleichskoeffizient für Feiertage eingefügt. Das zugrunde liegende Merkmal wird an Tagen nach einem Feiertag auf 1 gesetzt. Der dafür im ARIMA-Modell enthaltene Koeffizient soll den geringeren Verbrauch von Feiertagen ausgleichen. Dadurch verbessert sich das Modell auf den ersten Blick, da MAPE, Log-Likelihood und AIC etwas besser sind. Es handelt sich hierbei aber vermutlich um eine Überanpassung. Wie im Modell davor und im nächsten Modell erkennbar, führen beide Merkmale einzeln zu keiner Verbesserung.

#Skalierte Daten für Modellierung laden
df_scaled = pd.read_csv("../3-Data Preparation/data_scaled.csv", index_col=0, parse_dates=True)
df_scaled.index.freq = "D"

#Arbeitsfreie Tage erweitern
df_scaled["urlaubssaison"] = pd.read_csv("Daten/urlaubssaison_1.csv", index_col=0, parse_dates=True)
df_scaled.loc[df_scaled["urlaubssaison"] == 1, 'arbeitstag'] = 0

#Exogene Daten um Feiertag und Tag danach erweitern
df_scaled["feiertag"] = df_calendar["feiertag"].notna().astype(int)
df_scaled["feiertag_1"] = df_scaled["feiertag"].shift(1).fillna(0)

#Aufteilung in endogene und exogene Daten
exog_train = df_scaled[["arbeitstag","temperatur","feiertag","feiertag_1"]][:"2020-12-31"]
exog_test = df_scaled[["arbeitstag","temperatur","feiertag","feiertag_1"]]["2021-01-01":]

exog = exog_train.append(exog_test)

#Nicht-saisonale und saisonale Ordnung festlegen
order = (2, 0, 2) #p, d, q
seasonal_order = (2, 0, 2, 7) #P, D, Q, m

#Modell mit Trainingsdaten erstellen
train_model_config = SARIMAX(endog=endog_train, exog=exog_train, order=order, seasonal_order=seasonal_order)
train_model = train_model_config.fit()

#Modell mit allen Daten erstllen und Konfiguration/Koeffizienten von ersten Modell übernehmen (kein neues Training)
model_config = SARIMAX(endog=endog, exog=exog, order=order, seasonal_order=seasonal_order)
model = model_config.filter(train_model.params)

#Modell ausgeben
print(model.summary())
print()

#Vorhersage erzeugen
scaled_preds = model.predict(start="2021-01-1", end="2021-12-31", dynamic=False)

#Vorhersage invers-skalieren
scaler_target = joblib.load("../3-Data Preparation/scaler_endog.save")
preds = pd.DataFrame(data=scaler_target.inverse_transform(scaled_preds.values.reshape(-1, 1)), columns=["vorhergesagt"], index=pd.date_range('01/01/2021', periods=365, freq='D')).squeeze()

                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                         verbrauch   No. Observations:                 2557
Model:             SARIMAX(2, 0, 2)x(2, 0, 2, 7)   Log Likelihood                4670.613
Date:                           Thu, 03 Feb 2022   AIC                          -9315.227
Time:                                   22:26:40   BIC                          -9239.221
Sample:                               01-01-2015   HQIC                         -9287.665
                                    - 12-31-2021                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
arbeitstag     0.2092      0.002     95.737      0.000       0.205       0.213
temperatur    -0.0435      0.012     -3.485      0.000      -0.068      -0.019
feiertag      -0.1288      0.002    -52.764      0.000      -0.134      -0.124
feiertag_1    -0.0541      0.003    -21.477      0.000      -0.059      -0.049
ar.L1          0.8617      0.418      2.059      0.039       0.042       1.682
ar.L2          0.0435      0.378      0.115      0.908      -0.697       0.784
ma.L1         -0.0114      0.419     -0.027      0.978      -0.833       0.811
ma.L2         -0.0394      0.029     -1.365      0.172      -0.096       0.017
ar.S.L7        0.1424      0.059      2.412      0.016       0.027       0.258
ar.S.L14       0.8495      0.059     14.501      0.000       0.735       0.964
ma.S.L7       -0.0018      0.053     -0.035      0.972      -0.107       0.103
ma.S.L14      -0.7933      0.045    -17.824      0.000      -0.881      -0.706
sigma2         0.0015    2.1e-05     69.605      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                   1.10   Jarque-Bera (JB):              4724.07
Prob(Q):                              0.29   Prob(JB):                         0.00
Heteroskedasticity (H):               0.95   Skew:                            -0.24
Prob(H) (two-sided):                  0.47   Kurtosis:                         9.64
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

#Feiertage einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date, ymin=0, ymax=1, color="orange", linestyle="-.")

#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics(y_true, preds, True, False, True, False)

      Vorhersage
R2           1.0
MAE       3453.4
MSE   25427227.3
RMSE      5042.5
                
MAPE       2.0 %

Verbesserung durch Merkmal für vergangene Feiertage (nur 1)

In diesem Modell wird nur der Indikator für vergangene Feiertage (ohne den Indikator für Feiertage selbst) verwendet. Das Modell lässt sich dadurch allerdings nicht verbessern. Es ist also davon auszugehen, dass sich der Effekt der einzelnen Merkmale in Grenzen hält.

#Skalierte Daten für Modellierung laden
df_scaled = pd.read_csv("../3-Data Preparation/data_scaled.csv", index_col=0, parse_dates=True)
df_scaled.index.freq = "D"

#Arbeitsfreie Tage erweitern
df_scaled["urlaubssaison"] = pd.read_csv("Daten/urlaubssaison_1.csv", index_col=0, parse_dates=True)
df_scaled.loc[df_scaled["urlaubssaison"] == 1, 'arbeitstag'] = 0

#Exogene Daten um Feiertag und Tag danach erweitern
df_scaled["feiertag"] = df_calendar["feiertag"].notna().astype(int)
df_scaled["feiertag_1"] = df_scaled["feiertag"].shift(1).fillna(0)

#Aufteilung in endogene und exogene Daten
exog_train = df_scaled[["arbeitstag","temperatur","feiertag_1"]][:"2020-12-31"]
exog_test = df_scaled[["arbeitstag","temperatur","feiertag_1"]]["2021-01-01":]

exog = exog_train.append(exog_test)

#Nicht-saisonale und saisonale Ordnung festlegen
order = (2, 0, 2) #p, d, q
seasonal_order = (2, 0, 2, 7) #P, D, Q, m

#Modell mit Trainingsdaten erstellen
train_model_config = SARIMAX(endog=endog_train, exog=exog_train, order=order, seasonal_order=seasonal_order)
train_model = train_model_config.fit()

#Modell mit allen Daten erstllen und Konfiguration/Koeffizienten von ersten Modell übernehmen (kein neues Training)
model_config = SARIMAX(endog=endog, exog=exog, order=order, seasonal_order=seasonal_order)
model = model_config.filter(train_model.params)

#Modell ausgeben
print(model.summary())
print()

#Vorhersage erzeugen
scaled_preds = model.predict(start="2021-01-1", end="2021-12-31", dynamic=False)

#Vorhersage invers-skalieren
scaler_target = joblib.load("../3-Data Preparation/scaler_endog.save")
preds = pd.DataFrame(data=scaler_target.inverse_transform(scaled_preds.values.reshape(-1, 1)), columns=["vorhergesagt"], index=pd.date_range('01/01/2021', periods=365, freq='D')).squeeze()

                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                         verbrauch   No. Observations:                 2557
Model:             SARIMAX(2, 0, 2)x(2, 0, 2, 7)   Log Likelihood                4413.317
Date:                           Thu, 03 Feb 2022   AIC                          -8802.634
Time:                                   22:32:29   BIC                          -8732.475
Sample:                               01-01-2015   HQIC                         -8777.192
                                    - 12-31-2021                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
arbeitstag     0.3025      0.002    137.760      0.000       0.298       0.307
temperatur    -0.0618      0.014     -4.359      0.000      -0.090      -0.034
feiertag_1    -0.0309      0.002    -13.597      0.000      -0.035      -0.026
ar.L1          0.9433      1.698      0.556      0.579      -2.385       4.271
ar.L2         -0.0620      1.476     -0.042      0.967      -2.955       2.831
ma.L1         -0.0655      1.697     -0.039      0.969      -3.392       3.261
ma.L2         -0.0103      0.023     -0.452      0.651      -0.055       0.034
ar.S.L7        0.1604      0.296      0.541      0.588      -0.420       0.741
ar.S.L14       0.8300      0.295      2.816      0.005       0.252       1.408
ma.S.L7       -0.0440      0.304     -0.145      0.885      -0.640       0.552
ma.S.L14      -0.7130      0.265     -2.686      0.007      -1.233      -0.193
sigma2         0.0019   2.84e-05     66.073      0.000       0.002       0.002
===================================================================================
Ljung-Box (L1) (Q):                   1.31   Jarque-Bera (JB):              8890.99
Prob(Q):                              0.25   Prob(JB):                         0.00
Heteroskedasticity (H):               0.90   Skew:                             0.19
Prob(H) (two-sided):                  0.13   Kurtosis:                        12.13
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

#Feiertage einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date, ymin=0, ymax=1, color="orange", linestyle="-.")

#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics(y_true, preds, True, False, True, False)

      Vorhersage
R2           1.0
MAE       3559.7
MSE   26741693.1
RMSE      5171.2
                
MAPE       2.1 %

Verbesserung durch Merkmal für vergangene Feiertage (2)

Wenn zusätzlich zum Indikator für Feiertage zwei zusätzliche Koeffizienten für jeweils die nächsten beiden Tage zum Modell hinzugefügt werden, dann lässt sich das Modell nicht verbessern. Im Vergleich zu den vorherigen Modellen verschlechtern sich auch AIC und Log-Likelihood wieder.

#Skalierte Daten für Modellierung laden
df_scaled = pd.read_csv("../3-Data Preparation/data_scaled.csv", index_col=0, parse_dates=True)
df_scaled.index.freq = "D"

#Arbeitsfreie Tage erweitern
df_scaled["urlaubssaison"] = pd.read_csv("Daten/urlaubssaison_1.csv", index_col=0, parse_dates=True)
df_scaled.loc[df_scaled["urlaubssaison"] == 1, 'arbeitstag'] = 0

#Exogene Daten um Feiertag und Tage danach erweitern
df_scaled["feiertag"] = df_calendar["feiertag"].notna().astype(int)
df_scaled["feiertag_1"] = df_scaled["feiertag"].shift(1).fillna(0)
df_scaled["feiertag_2"] = df_scaled["feiertag"].shift(2).fillna(0)

#Aufteilung in endogene und exogene Daten
exog_train = df_scaled[["arbeitstag","temperatur","feiertag","feiertag_1","feiertag_2"]][:"2020-12-31"]
exog_test = df_scaled[["arbeitstag","temperatur","feiertag","feiertag_1","feiertag_2"]]["2021-01-01":]

exog = exog_train.append(exog_test)

#Nicht-saisonale und saisonale Ordnung festlegen
order = (2, 0, 2) #p, d, q
seasonal_order = (2, 0, 2, 7) #P, D, Q, m

#Modell mit Trainingsdaten erstellen
train_model_config = SARIMAX(endog=endog_train, exog=exog_train, order=order, seasonal_order=seasonal_order)
train_model = train_model_config.fit()

#Modell mit allen Daten erstllen und Konfiguration/Koeffizienten von ersten Modell übernehmen (kein neues Training)
model_config = SARIMAX(endog=endog, exog=exog, order=order, seasonal_order=seasonal_order)
model = model_config.filter(train_model.params)

#Modell ausgeben
print(model.summary())
print()

#Vorhersage erzeugen
scaled_preds = model.predict(start="2021-01-1", end="2021-12-31", dynamic=False)

#Vorhersage invers-skalieren
scaler_target = joblib.load("../3-Data Preparation/scaler_endog.save")
preds = pd.DataFrame(data=scaler_target.inverse_transform(scaled_preds.values.reshape(-1, 1)), columns=["vorhergesagt"], index=pd.date_range('01/01/2021', periods=365, freq='D')).squeeze()

                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                         verbrauch   No. Observations:                 2557
Model:             SARIMAX(2, 0, 2)x(2, 0, 2, 7)   Log Likelihood                4621.884
Date:                           Thu, 03 Feb 2022   AIC                          -9215.769
Time:                                   22:26:51   BIC                          -9133.917
Sample:                               01-01-2015   HQIC                         -9186.087
                                    - 12-31-2021                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
arbeitstag     0.2292      0.002    101.048      0.000       0.225       0.234
temperatur     0.0080      0.013      0.635      0.525      -0.017       0.033
feiertag      -0.1051      0.003    -37.919      0.000      -0.111      -0.100
feiertag_1    -0.0486      0.003    -15.711      0.000      -0.055      -0.042
feiertag_2     0.0223      0.003      7.084      0.000       0.016       0.028
ar.L1          0.7617      1.382      0.551      0.582      -1.947       3.471
ar.L2          0.1609      1.289      0.125      0.901      -2.365       2.686
ma.L1          0.1551      1.382      0.112      0.911      -2.553       2.863
ma.L2         -0.0143      0.025     -0.578      0.563      -0.063       0.034
ar.S.L7        0.2445      0.178      1.373      0.170      -0.105       0.594
ar.S.L14       0.7339      0.176      4.180      0.000       0.390       1.078
ma.S.L7       -0.0722      0.172     -0.420      0.675      -0.409       0.265
ma.S.L14      -0.6515      0.140     -4.638      0.000      -0.927      -0.376
sigma2         0.0015   2.26e-05     66.515      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                   2.55   Jarque-Bera (JB):              4108.44
Prob(Q):                              0.11   Prob(JB):                         0.00
Heteroskedasticity (H):               0.94   Skew:                            -0.08
Prob(H) (two-sided):                  0.41   Kurtosis:                         9.21
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

#Feiertage einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date, ymin=0, ymax=1, color="orange", linestyle="-.")

#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics(y_true, preds, True, False, True, False)

      Vorhersage
R2           1.0
MAE       3519.0
MSE   25980606.0
RMSE      5097.1
                
MAPE       2.0 %

Verbesserung durch Merkmal für vergangene Feiertage (7)

Das Modell lässt sich auch mit einem zusätzlichen Ausgleichskoeffizienten für den siebten Tag nach einem Feiertag nicht mehr verbessern.

#Skalierte Daten für Modellierung laden
df_scaled = pd.read_csv("../3-Data Preparation/data_scaled.csv", index_col=0, parse_dates=True)
df_scaled.index.freq = "D"

#Arbeitsfreie Tage erweitern
df_scaled["urlaubssaison"] = pd.read_csv("Daten/urlaubssaison_1.csv", index_col=0, parse_dates=True)
df_scaled.loc[df_scaled["urlaubssaison"] == 1, 'arbeitstag'] = 0

#Exogene Daten um Feiertag und Tage danach erweitern
df_scaled["feiertag"] = df_calendar["feiertag"].notna().astype(int)
df_scaled["feiertag_1"] = df_scaled["feiertag"].shift(1).fillna(0)
df_scaled["feiertag_2"] = df_scaled["feiertag"].shift(2).fillna(0)
df_scaled["feiertag_7"] = df_scaled["feiertag"].shift(7).fillna(0)

#Aufteilung in endogene und exogene Daten
exog_train = df_scaled[["arbeitstag","temperatur","feiertag","feiertag_1","feiertag_2","feiertag_7"]][:"2020-12-31"]
exog_test = df_scaled[["arbeitstag","temperatur","feiertag","feiertag_1","feiertag_2","feiertag_7"]]["2021-01-01":]

exog = exog_train.append(exog_test)

#Nicht-saisonale und saisonale Ordnung festlegen
order = (2, 0, 2) #p, d, q
seasonal_order = (2, 0, 2, 7) #P, D, Q, m

#Modell mit Trainingsdaten erstellen
train_model_config = SARIMAX(endog=endog_train, exog=exog_train, order=order, seasonal_order=seasonal_order)
train_model = train_model_config.fit()

#Modell mit allen Daten erstllen und Konfiguration/Koeffizienten von ersten Modell übernehmen (kein neues Training)
model_config = SARIMAX(endog=endog, exog=exog, order=order, seasonal_order=seasonal_order)
model = model_config.filter(train_model.params)

#Modell ausgeben
print(model.summary())
print()

#Vorhersage erzeugen
scaled_preds = model.predict(start="2021-01-1", end="2021-12-31", dynamic=False)

#Vorhersage invers-skalieren
scaler_target = joblib.load("../3-Data Preparation/scaler_endog.save")
preds = pd.DataFrame(data=scaler_target.inverse_transform(scaled_preds.values.reshape(-1, 1)), columns=["vorhergesagt"], index=pd.date_range('01/01/2021', periods=365, freq='D')).squeeze()

                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                         verbrauch   No. Observations:                 2557
Model:             SARIMAX(2, 0, 2)x(2, 0, 2, 7)   Log Likelihood                4636.493
Date:                           Thu, 03 Feb 2022   AIC                          -9242.986
Time:                                   22:27:03   BIC                          -9155.287
Sample:                               01-01-2015   HQIC                         -9211.184
                                    - 12-31-2021                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
arbeitstag     0.2156      0.002     94.311      0.000       0.211       0.220
temperatur    -0.0404      0.013     -3.088      0.002      -0.066      -0.015
feiertag      -0.1239      0.003    -44.387      0.000      -0.129      -0.118
feiertag_1    -0.0567      0.003    -17.087      0.000      -0.063      -0.050
feiertag_2     0.0085      0.004      2.366      0.018       0.001       0.016
feiertag_7     0.0019      0.004      0.495      0.621      -0.006       0.009
ar.L1          0.8002      0.918      0.872      0.383      -0.998       2.599
ar.L2          0.1206      0.851      0.142      0.887      -1.547       1.788
ma.L1          0.1639      0.917      0.179      0.858      -1.634       1.962
ma.L2         -0.0139      0.040     -0.346      0.729      -0.093       0.065
ar.S.L7        0.2302      0.186      1.238      0.216      -0.134       0.595
ar.S.L14       0.7508      0.184      4.091      0.000       0.391       1.111
ma.S.L7       -0.0659      0.180     -0.367      0.714      -0.418       0.286
ma.S.L14      -0.6681      0.148     -4.517      0.000      -0.958      -0.378
sigma2         0.0015   2.34e-05     65.587      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                  18.60   Jarque-Bera (JB):              4097.87
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               0.93   Skew:                            -0.08
Prob(H) (two-sided):                  0.26   Kurtosis:                         9.20
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

#Feiertage einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date, ymin=0, ymax=1, color="orange", linestyle="-.")

#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics(y_true, preds, True, False, True, False)

      Vorhersage
R2           1.0
MAE       3471.5
MSE   25663303.2
RMSE      5065.9
                
MAPE       2.0 %

Verbesserung durch Merkmal für vergangene Feiertage (14)

Letztlich wird noch ein Ausgleichskoeffizient für den vierzehnten Tag nach einem Feiertag hinzugefügt. Auch hierdurch ist keine weitere Verbesserung zu erzielen.

#Skalierte Daten für Modellierung laden
df_scaled = pd.read_csv("../3-Data Preparation/data_scaled.csv", index_col=0, parse_dates=True)
df_scaled.index.freq = "D"

#Arbeitsfreie Tage erweitern
df_scaled["urlaubssaison"] = pd.read_csv("Daten/urlaubssaison_1.csv", index_col=0, parse_dates=True)
df_scaled.loc[df_scaled["urlaubssaison"] == 1, 'arbeitstag'] = 0

#Exogene Daten um Feiertag und Tage danach erweitern
df_scaled["feiertag"] = df_calendar["feiertag"].notna().astype(int)
df_scaled["feiertag_1"] = df_scaled["feiertag"].shift(1).fillna(0)
df_scaled["feiertag_2"] = df_scaled["feiertag"].shift(2).fillna(0)
df_scaled["feiertag_7"] = df_scaled["feiertag"].shift(7).fillna(0)
df_scaled["feiertag_14"] = df_scaled["feiertag"].shift(14).fillna(0)

#Aufteilung in endogene und exogene Daten
exog_train = df_scaled[["arbeitstag","temperatur","feiertag","feiertag_1","feiertag_2","feiertag_7","feiertag_14"]][:"2020-12-31"]
exog_test = df_scaled[["arbeitstag","temperatur","feiertag","feiertag_1","feiertag_2","feiertag_7","feiertag_14"]]["2021-01-01":]

exog = exog_train.append(exog_test)

#Nicht-saisonale und saisonale Ordnung festlegen
order = (2, 0, 2) #p, d, q
seasonal_order = (2, 0, 2, 7) #P, D, Q, m

#Modell mit Trainingsdaten erstellen
train_model_config = SARIMAX(endog=endog_train, exog=exog_train, order=order, seasonal_order=seasonal_order)
train_model = train_model_config.fit()

#Modell mit allen Daten erstllen und Konfiguration/Koeffizienten von ersten Modell übernehmen (kein neues Training)
model_config = SARIMAX(endog=endog, exog=exog, order=order, seasonal_order=seasonal_order)
model = model_config.filter(train_model.params)

#Modell ausgeben
print(model.summary())
print()

#Vorhersage erzeugen
scaled_preds = model.predict(start="2021-01-1", end="2021-12-31", dynamic=False)

#Vorhersage invers-skalieren
scaler_target = joblib.load("../3-Data Preparation/scaler_endog.save")
preds = pd.DataFrame(data=scaler_target.inverse_transform(scaled_preds.values.reshape(-1, 1)), columns=["vorhergesagt"], index=pd.date_range('01/01/2021', periods=365, freq='D')).squeeze()

                                     SARIMAX Results                                     
=========================================================================================
Dep. Variable:                         verbrauch   No. Observations:                 2557
Model:             SARIMAX(2, 0, 2)x(2, 0, 2, 7)   Log Likelihood                4652.339
Date:                           Thu, 03 Feb 2022   AIC                          -9272.679
Time:                                   22:27:16   BIC                          -9179.133
Sample:                               01-01-2015   HQIC                         -9238.756
                                    - 12-31-2021                                         
Covariance Type:                             opg                                         
===============================================================================
                  coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------------------------------------------
arbeitstag      0.2090      0.002     90.261      0.000       0.204       0.214
temperatur     -0.0433      0.013     -3.304      0.001      -0.069      -0.018
feiertag       -0.1354      0.003    -47.666      0.000      -0.141      -0.130
feiertag_1     -0.0551      0.003    -17.234      0.000      -0.061      -0.049
feiertag_2      0.0077      0.004      2.141      0.032       0.001       0.015
feiertag_7      0.0009      0.004      0.217      0.828      -0.007       0.009
feiertag_14    -0.0035      0.004     -0.850      0.395      -0.011       0.005
ar.L1           0.8298      0.337      2.465      0.014       0.170       1.490
ar.L2           0.0858      0.308      0.278      0.781      -0.518       0.690
ma.L1           0.0154      0.338      0.046      0.964      -0.647       0.678
ma.L2          -0.0594      0.029     -2.083      0.037      -0.115      -0.004
ar.S.L7         0.1832      0.052      3.495      0.000       0.080       0.286
ar.S.L14        0.8044      0.052     15.532      0.000       0.703       0.906
ma.S.L7        -0.0188      0.046     -0.404      0.686      -0.110       0.072
ma.S.L14       -0.7693      0.038    -20.375      0.000      -0.843      -0.695
sigma2          0.0015   2.35e-05     65.260      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                   0.36   Jarque-Bera (JB):              4510.03
Prob(Q):                              0.55   Prob(JB):                         0.00
Heteroskedasticity (H):               0.95   Skew:                            -0.23
Prob(H) (two-sided):                  0.41   Kurtosis:                         9.49
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

#Feiertage einzeichnen
for date in df_calendar["2021-01-01":][df_calendar["feiertag"].isna() == False].index:

    plt.axvline(x=date, ymin=0, ymax=1, color="orange", linestyle="-.")

#Vorhersagen mit echten Werten vergleichen
functions.custom_metrics(y_true, preds, True, False, True, False)

      Vorhersage
R2           1.0
MAE       3496.4
MSE   26123208.8
RMSE      5111.1
                
MAPE       2.0 %

Ergebnis

Stellenweise lässt sich das ARIMA-Modell durch die Ausgleichskoeffizienten verbessern, die Verbesserungen sind allerdings gering ausgeprägt und führen teilweise auch zu Überanpassungen. Ein wesentliches Problem ist, dass es nur etwa 30 Feiertage mit dem entsprechenden Problem in den Daten gibt. Bei insgesamt 2.557 Datensätzen sind daher sehr wenig Beispieldaten vorhanden, anhand derer der Algorithmus Daten generalisieren kann. Warum auf Oversampling verzichtet wird, ist in der Ausarbeitung genauer dargestellt. Vereinfacht gesagt lässt sich die Zeitreihe nicht ohne weiteres durch zusätzliche Beobachtungen erweitern, ohne das Wesen und den Verlauf der Zeitreihe zu ändern (Saisonalitäten, Autokorrelationen, Zusammenhänge mit anderen Merkmalen etc.). Unabhängig davon tritt das Problem nicht bei allen Feiertagen auf, folgt beispielsweise ein arbeitsfreier Tag auf den Feiertag ist das Problem sehr viel geringer, hier wären dann Ausgleichskoeffizienten für die Ausgleichskoeffizienten nötig. Der Grundgedankt, alle möglichen Sonderfälle durch eigene Koeffizienten abzudecken, widerspricht allerdings dem Grundgedanken des maschinellen Lernens: Der Generalisierung von Informationen aus vorhandenen Daten. Weiterhin ist zu erwähnen, dass das Modell sowohl die Baseline als auch das Erfolgskriterium bereits (mit Abstand) schlagen kann. Aus diesen Gründen wird vom Einsatz der hier beschriebenen Merkmale abgesehen.